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Abstract 

In this work we apply the Landau-Lifshitz pseudotensor flux formal¬ 
ism to the calculation of the total mass and the total angular momentum 
during the evolution of a binary black hole system. We also compare its 
performance with the traditional integrations for the global quantities. 

It shows that the advantage of the pseudotensor flux formalism is the 
smoothness of the numerical value of the global quantities, especially of 
the total angular momentum. Although the convergence behavior of the 
global quantities with the pseudotensor flux method is only comparable 
with the ones with the traditional method, the smoothness of its numer¬ 
ical value allows using a larger radius for surface integration to obtain 
more accurate result. 

This work is dedicated to the General Relativity special issue. 

PACS numbers: 04.25.D-, 04.25.dg, 04.25.Nx, 04.70.-s 


1 Introduction 

Since the breakthrough by Pretorius [I] in 2005, the inspiral, merger, and ring- 
down of binary compact objects (including binary black holes, i.e., BBHs, binary 
neutron stars, black hole-neutron star binary) has been successfully simulated 
to high accuracy. The emphasis in the field is now turned to extracting as- 
trophysical information from these simulations. Therefore, it will be useful in 
validation to measure one physical quantity with different methods. 

The pseudotensor formalisms [5] and the quasilocal quantities mm arise 
from viewing general relativity as a nonlinear field theory in a fixed background 
reference, especially in a flat auxiliary spacetime. These formalisms have been 
used to explore the nonlinear dynamics of spacetime, for example, dynamical 
horizons , the distribution and flow of linear momentum in stronly nonlinearly 
curved spacetimes [5], and black hole spin measurement [7]. The result shows 
the usefulness of these formalisms and also shed a light on their possibly broader 
applications as analytical tools in various numerical simulations. 
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In this work we would like to try an alternative method in calculating the to¬ 
tal mass and the total angular momentum with the Landau-Lifshitz pseudoten¬ 
sor formalism. The motivation comes from the imperfection of the commonly 
used method, i.e., Eqs. (IMl) and (I55|) . in calculating the Arnowitt-Deser-Misner 
(ADM) mass and the angular momentum. The calculation of these two global 
quantities basically executes an integration over a two-sphere on the spatial do¬ 
main. Due to the limit of the grid resolution, the numerical values of these two 
quantities, especially the one of the angular momentum, fluctuate around the 
average values. The numerical fluctuation makes it difficult to tell the numerical 
value accurately. Different from the ADM mass and the angular momentum, 
the calculation of the momentum flux from the Landau-Lifshitz pseudotensor 
formalism includes not only an integration over a two-sphere on the spatial do¬ 
main, but also an integration over the time domain. We expect that this method 
with the extra integration will give smoother global quantity curves with respect 
to time, and at least as accurate as the one for the ADM mass and the angular 
momentum. 

The rest of this work is organized as follows: In the next section, we give a 
description of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation 
commonly used in numerical relativity. We then describe the methods used 
in this work for the calculation of the global quantities, i.e., the total mass 
and the angular momentum, in Sec. [3] In this section, besides the usual ADM 
and the angular momentum calculation, it gives the integral formulas of these 
global quantities with the Landau-Lifshtiz pseudotensor. We then report the 
numerical comparison between these two different methods in the simulation of 
BBH with spins in Sec. SI And the discussion and summary will be presented 
in the Sec. [3 Throughout the paper, geometric units with G = c = 1 are used. 
Einstein summation rule is adopted unless stated explicitly. 


2 The BSSN Formulation 

The metric in the ADM form is 

ds^ = —a^dt^ + (dx* -I- /3Mt)(dx'^ -I- jd^dt), (1) 

wherein a is the lapse function, /3* is the shift vector, and is the spatial 
three-metric. Throughout this paper, Latin indices are spatial indices and run 
from 1 to 3, whereas Greek indices are space-time indices and run from 0 to 3. 

Einstein’s equations can then be decomposed into the Hamiltonian constraint 
T-L and the momentum constraints Aii 

n = R- +k'^ = ( 2 ) 

M^ = = 0, (3) 

and the evolution equations 

= —‘2aKij, 
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( 4 ) 


( 5 ) 



—'S/iS/jOi + a^Rij — 2Ki£K^j + KKij). 


Here we have assumed vacuum Tap = Q and have used 


d 

dt 


dt 




( 6 ) 


where £^ is the Lie derivative with respect to Vi is the covariant derivative 
associated with jij , Rij is the three-dimensional Ricci tensor 

Rij — 2^ ^klpj ^ij,ki) “t“ L iTmkj £mij£ fcj (7) 

where ^ 

+ l^k,j - Ijk.i)- ( 8 ) 

And R is its trace R = Rij. 

In the BSSN formalism [S], the above ADM equations are rewritten by in¬ 
troducing the conformally related metric 7 ij 

( 9 ) 


with the conformal exponent cj) chosen so that the determinant 7 of ^ij is unity 


= 7'/^ 


( 10 ) 


where 7 is the determinant of jij. The traceless part of the extrinsic curvature 
Kij, defined by 

Aij = Al(ij) = Kij — —'jijK, (11) 

where Kij with two indices between () is to take the symmetric and traceless part 
of Kij, and K = j^^Kij is the trace of the extrinsic curvature, is conformally 
decomposed according to 

A, = ( 12 ) 

The conformal connection functions T*, initially defined as 

r = A''r% = -f\p ( 13 ) 

are regarded as independent variables in this formulation. 

The evolution equations of BSSN formulation can be written as 


d _ - 

^ 7 ij = —2,aAij, 

= a + ^K^^ - V^a, 


(14) 

(15) 

(16) 
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(17) 


—Ay = a{KAij — 2AikA^j) + e — V(iVj)a), 

5tf* = 2a (f]kA^'^ - - 2A*^aj 

+ + \i^^p\,k. (18) 

The Ricci tensor i?y can be written as a sum of two pieces 

%=%+4, (19) 


where is given by 

Rf^ = -2V,Vj(j) - + 4V,0Vj(/) - 47yV'=(/)Vfc(/), (20) 

where is the covariant derivative with respect to 7 y , while, with the help of 


the r*, Rij can be expressed as 

Rij = ,j) + f^f (y)fe + 2f^^(jf iVklj ■ (21) 

The new variables are tensor densities, so that their Lie derivatives are 

£^K = P^K^k, ( 22 ) 

£^P = P'^P,k + \p\k, (23) 

£pli3 = P^lij,k + 2pk(iP^,3) - (24) 

£0-^ij ~ + 2Afc(j/3^y-) — —AijP^^k- (25) 

The Hamiltonian and momentum constraints ([2]) and © can be rewritten as 

n = e-‘^'^iR - 8V^p - 8VVV,(/)) + - iyi^ = 0, (26) 

M^ = VjAp + ep^jAp - = 0, (27) 


where R = 7 *-^ i?y. 


3 Calculation of Global Quantities 

3.1 ADM mass and angular momentum 

The ADM mass is defined in terms of a surface integral at spatial infinity. In nu¬ 
merical simulations, this integral can be approximated by an integral evaluated 
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on a surface near the outer boundaries of the grid. In Cartesian coordinates, 
the ADM mass is defined by a surface integral at spatial infinity [3110] 




(28) 


where d^E^ = {l/2)y/^eijkdx ^is the surface element and e^-fe the Levi-Civita 
alternating symbol. We now perform a conformal decomposition 

(29) 

where ijj = e^. Assuming the asymptotic behavior 

t/j ^ 1 + 0{-) when r —>■ oo, (30) 

r 


and 


we can rewrite 


7 b ^ + O(-) when r —>■ oo, 

as 


(31) 


'll: ['(/’^(ymn.i - 7in,m) + 4V'^('l/'.i7mn - V'.mbjn)] d^Ei 


1 

IOtt 


J OO 


_L / (r - - 8VV)d^S» = 7 ^ (f* - 8V*e'^)d2S,. 


IGtt 


IOtt /oo 


(32) 


Here the conformal surface element is defined as d^E^ = (l/2)eyfcdx'’dx^ since 
7 = 1, we use the abbreviations f* = and = 0, and 

is the three-covariant derivative with respect to the metric . 

We define the angular momentum J* as (compare [TTl IT^ l 

/ = / x^A^fcd^E, = / x^e®‘^i'fcd"E,, (33) 

Stt /oo Stt " /oo 

where the indices of 6^^ are raised and lowered with the flat metric 5ij, d^E^ = 
e^’^d^S,, and A^j = A^j. 

Therefore, the surface integrals of the ADM mass and the angular momen¬ 
tum (in vacuum) are respectively m- 

M =— / (r - 8W^e*)d^ti, (34) 

Idtt Jan 

J^ (35) 

Stt /an 

These two global quantities are useful tools for the system diagnostics to validate 
the calculations. 
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3.2 Pseudotensor and momentum flux 

In this section, we briefly review the Landau-Lifshitz formulation of gravity 
and the statement of four-momentum conservation within this theory. The 
Landau-Lifshitz formulation has been described in HE] to reformulate general 
relativity as a nonlinear field theory in flat spacetime. Here we follow closely 
to the content in H. In this formalism, an arbitrary asymptotically Lorentz 
coordinate is firstly built on a given curved (but asymptotically-flat) spacetime. 
Then the coordinate is used to map the curved (i.e. physical) spacetime onto 
an auxiliary flat spacetime by enforcing that the coordinate on this spacetime 
are globally Lorentz. The auxiliary flat metric takes the Minkowski form, rj^i, = 
diag(-l, 1,1,1). 

Gravity is described, in this formulation, by the physical metric density 

(36) 

where g is the determinant of the covariant components of the physical metric, 

and are the contravariant components of the physical metric. In terms of 

the superpotential 

(37) 

the Einstein field equations take the field-theory-in-flat-spacetime form 

= IGttt'^^ (38) 

Here = {—g){T^'^ + t^) is the total effective stress-energy tensor, indices 
after the comma denote partial derivatives (covariant derivatives with respect 
to the flat auxiliary metric), and the Landau-Lifshitz pseudotensor (actually 
a real tensor in the auxiliary flat spacetime) is given by 

I67r{-g)t^l = + ^g^’^gx.Q^^pQ^'^,r 

+ ^i2g^^g’^^ - g^'^/")(25.pg., - 5p«ff.^)0"^A0^",a. (39) 

And the Landau-Lifshitz pseudotensor can be expressed in term of the 4-metric 
and the 4-connection as 

i67rt^L - L^) + 2r’^^''L^ - (r^ - L^)(r'^ - l'') + 

- Fxa^r^^'' + g>^''{L,L'^ - 2L,r^ + FxapF^^P), 

(40) 

where F^ = g^'^F^^xa, = F’^^cr- With the relation equations in Appendix lAl 

the equation can be re-expressed easily with the 3-1-1 quantities. By virtue of the 
symmetries of the superpotential (which are the same as those of the Riemann 
tensor), the field equations in the form (l38l) imply the differential conservation 
law for four-momentum 

= 0, (41) 
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which is equivalent to = 0 (where the semicolon denotes a covariant 

derivative with respect to the physical metric). 

It is shown in w that the total four-momentum of any isolated system as 
measured gravitationally in the asymptotically flat region far from the system 
is 

(42) 

Jv 

Thus 

= ^J = J = - J = - ^ , (43) 

where d^S^ = (l/2)eyfcda;'’da;^ is the surface-area element (defined by using the 
flat auxiliary metric), and the integral is over an arbitrarily large closed surface 
S surrounding the system, and Eq. (1411) is used. Therefore, the total momentum 
flux across the 2-surface within [^ 1 ,^ 2 ] is 

^Ptot = Ptotih) - Ptotiti) = - [ <f (44) 

Jti Js 

With = M, this leads to 

M{t) = M(0) - [ (f (45) 

Jo Js 

where M(0) can be obtained by using Eq. (IMl) at t = 0. 

The total angular momentum of any isolated system as measured gravita¬ 
tionally in the asymptotically flat region far from the system is 


JZ = 2 / X^>^T''^°d^X. 


(46) 


Thus 

dJZ ^ 2 


xK’'^°d^x = 2 J {x^>^Z°)^od^x = 2 J +x'^>^Z°.p)d^x 

= 2 [ - {x^>^Z^)^j]dZ = - (f {x^t'^^ - x'^T^^)d^tj, (47) 

Jv ’ Js 

Therefore, the total angular momentum flux across the 2-surface S within [ti, t 2 ] 
is 

^ JZtih) - JZtiti) = - r <{ {x^^r’'^ - x’'T>^nd%dt, (48) 

Jti J s 


and for 


J^(t) = J^(0) — f (f {xZ^ — yZ^)d^'E,jdt. (49) 

Jo Js 

where T^(0) can be obtained by using Eq. (l3^ at t = 0. 

We will use Eqs. and (l4^ to calculate the mass and the angular mo¬ 
mentum, and compare them with the result from Eqs. (I34p and (1351) . 
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4 Numerical Result 


The AMSS-NCKU code with the standard moving box style mesh refinement 
[Tl[T5l|T6] is used in this work. We used 10 mesh levels, and the finest 3 levels 
are movable in evolving the binary black holes (BBHs). In each fixed level, we 
used one box with 128 x 128 x 64 grids with assumed equatorial symmetry. The 
outermost physical boundary is 512M and this makes the finest resolution to be 
h = M/64. For the movable levels, two boxes with 64 x 64 x 32 grids are used 
to cover each black hole. In time direction, the Berger-Oliger numerical scheme 
is adopted for the levels higher than four. 

The moving puncture gauge condition 


dta = — 2aK, 

(50) 


(51) 

dtB^ = dtr - r]B^ + P^B^j - 

(52) 


is used and has been shown to give good behavior for the black hole simulations 
in [14] . In this paper we use f] = 2M with M being the ADM mass of the given 
configuration. 

In this section we apply the analysis tools described in the above section to 
the inspiralling binary black hole systems. We present two cases in this paper. 
One corresponds to the spinless binary black holes initially. The another one 
corresponds to two fast-spinning black holes initially. The two individual black 
holes in the binary are identical in the both cases. In the fast-spinning case, 
the spin parameter for each black hole is a = 0.9. And the spin is aligned 
to the orbital angular momentum. For the detailed description of the initial 
data construction, the grid setting for the numerical evolution, and the involved 
numerical tricks, we refer our reader to m- 

In Fig.[T]we compare the binary’s mass and its angular momentum calculated 
with the traditional integrations, i.e., Eqs. (IM)) and (l35ll . and the pseudotensor 
flux integrations, i.e, Eqs. (gHI) and (H^ . We show the results for three differ¬ 
ent extraction radii r = 50, 80, and 120 respectively. In the both BBH cases, 
the physical quantities calculated with the traditional integration allow larger 
fluctuations which come from the numerical error. For the convenience of com¬ 
parison, we smooth the traditional data by averaging within each time range 
5M. In the figure, we denote the data as “ADM after smooth”. And the data 
marked with “ADM” corresponds to the raw data. The data after smoothing 
becomes much smoother. However, by the comparison with the quantities from 
the pseudotensor flux integrations, the smoothed data still fluctuates more. In 
the plot of mass, such fluctuation appears after the junk radiation reaches the 
extraction sphere. We consider such fluctuation as the gauge adjustment re¬ 
sulted from the junk radiation. The numerical error could also come from the 
reflection of the junk radiation via the mesh refinement boundary, and thus 
contribute to the fluctuation. As we can see from Fig. (T] for the traditional 
integration, the angular momentum is even more sensitive to these factors. So 
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r = 50 r = 50 




r = 80 




r= 120 


r= 120 




Figure 1: Comparisons of the physical quantities calculated with the traditional 
integration and the pseudotensor flux integrations for the spinless BBH case. 
The left column corresponds to the masses M measured at the radii r = 50, 80, 
and 120 respectively. The right column corresponds to the angular momenta Jz 
of the spacetime measured at the same radii as in the left column. 
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r=50 



r = 80 



r= 120 



Figure 2: Same as in Fig. [U except that these plots are for the BBH case with 
spin a = 0.9. 


r = 50 




T[M] 


r= 120 
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Figure 3: Comparison of the effect of the finite extraction radius on the value 
of mass calculated with the traditional integration and the pseudotensor flux 
integration for the spinless BBH case. The data for the ADM masses here have 
been smoothed as explained in Fig. [TJ 


even after the junk radiation passes away, the angular momentum still fluctuates 
mildly due to the numerical error from the mesh refinement boundary reflection. 
Interestingly, the quantities calculated with the pseudotensor flux integrations 
seem immune to these factors. Figure gives the similar result, but for fast¬ 
spinning BBH case. For the fast-spinning BBH case, the gauge dynamics is 
more complicated. So we can see the fluctuation of the traditional integration 
is more drastic than the spinless case. The result from the pseudotensor flux 
calculation still works smoothly in this extreme configuration. Except those 
fluctuation in the traditional integration, the results of these two analysis tools 
are consistent to each other in both Fig. [T] and [2] 

The global quantities are formally defined at infinity. However, we can only 
calculate them at some finite radius in practice. This may cause some ambiguity. 
In principle, the sequence corresponding to different extraction radii should con¬ 
verge to the quantities defined at infinity. To be a good analysis tool, we expect 
that the method gives a fast convergence. In Fig. [3] we compare the convergence 
behavior of the mass integral with respect to different extraction radii with the 
traditional integration and the pseudotensor flux integration for the spinless 
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Figure 4: Same as in Fig. |3l except that these plots are for the BBH case with 
spin a = 0.9. 


BBFl case. During the junk radiation period, the convergence of the pseudoten¬ 
sor flux method is roughly two times better than the traditional integration. But 
during the merger part, the convergence of the traditional integration method is 
two times better than the pseudotensor flux method. Considering that the junk 
radiation is unphysical, we conclude that the traditional integration method is 
a better analysis tool in this aspect. In Fig. 01 we did the same investigation 
for fast-spinning BBH case. For the junk radiation part the same result can 
be seen as the spinless case. For the merger part, the fast-spinning BBH con¬ 
figuration introduces some challenge to the numerical evolution as explained in 
m- So as ones expect, the convergence behaviors for both analysis methods 
are equally bad, although they are consistent with each other. As to the angular 
momentum, the result from the traditional integration fluctuate so much that 
it does not make sense to compare it with the one from the pseudotensor flux 
integration. 


5 Summary 

In this work we apply the Landau-Lifshitz pseudotensor flux formalism as an 
alternative method in calculating the total mass and the total angular momen- 
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turn during the evolutions of a binary black hole system. We also compare its 
performance with the traditional integrations for the global quantities. Due to 
the gauge choice employing a flat spacetime background in the Landau-Lifshitz 
pseudotensor formalism, it is not expected that the result from this method will 
be accurate enough for the radius for integration is not far from the singularity. 
However, we find that the overall result with the method is consistent with the 
one with the traditional integration. 

The advantage of the pseudotensor flux formalism is the smoothness of the 
global quantities, especially of the total angular momentum. It has been plagued 
for a long time with the fluctuation and inaccuracy of the numerical value of the 
total angular momentum calculated with the traditional integration, especially 
when the grid resolution is usually low for the radius of the surface integration 
is large. It shows in this work that this problem can be solved with the pseu¬ 
dotensor flux method. The reason mainly comes from the integrations along 
the time domain in Eqs. (HSI) and (HU). Therefore, although the convergence 
behavior of the global quantities with the pseudotensor flux method is only 
comparable with the ones with the traditional method, the smoothness of its 
numerical value allows using a larger radius for surface integration to obtain 
more accurate result. 

As showed in [TB] and HZ], the total angular momentum calculated with the 
traditional method usually decays after the merger in the fast-spinning BBH 
cases. In our BBH simulations, it seems that the total angular momentum with 
the pseudotensor flux method conserves much better than the one with the 
traditional method. However, it might need a further detailed investigation to 
confirm this point. 

This work shows that the pseudotensors (and the quasi-local quantities) 
could be very useful analysis tools in numerical relativity. Therefore, we plan 
to study the usefulness of different pseudotensors/quasi-local quantities, and 
also the advantage of different spacetime background, e.g, the Schwarzschild 
spacetime or the Kerr spacetime, in numerical relativity in the future. 
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A Spacetime 4-connection in 3-1-1 Expression 

The 4-inetric can be constructed out of the 3-metric 7 ^- and the lapse a and 
shift functions /3® as 


5oo 

90k 


P^PI- 

Pk 

9]0 

9jk 



Ijk 

0 

0 

gOk 


1 

pk 


a2 

a2 

5 ^° 

gjk 



' 2 


(53) 


(54) 


where Pi = jijP^ ■ 

From Appendix B of [TO] we can obtain the following expressions for the 
4-connection in terms of 3-1-1 quantities 


r° - = -—K 


a 


1 


= -(V,a - K^^pn = V, Ina + r\^P'^, 
a 

r%o = -{dta + p-^Vma - K^nP^^P^ = Ot \n a + 
a 


Omi 


pi _ “pz 

^ jk — ^ jk 


P^ 


jk. 


P' 


'i pO 


r\, = Vj/3* - aK^, + - Vjtt) = V,/3* - aK^j - P^F 

r\o = dtP^ + P^V^P^ + a{V^a - 2iFV/3™) 


Oj, 


(55) 

(56) 

(57) 

(58) 

(59) 


+ - dta - /3™V^a) 

a 

= dtP^ + p^r\m + (a' 7 *™ + P^P^)r\m - P^r%o, 


(60) 


where Vi is the covariant derivative associated with the 3-metric 7 ij, and the 
corresponding 3-connection F^^fe. For 


Fo = 5tlna + V^/3'"-aiF, 

F, = 5ana + r’”^„ 

F° = ^{p^d^a-dta-a^K), 


(61) 

(62) 

(63) 


r = r- i(d^a - P^K) - \{dtp^ - P^d^P^) + ^{dta - P^d^a), (64) 


where F® = ^^^T^jk- 
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